Vortex Rings and Lieb Modes in a Cylindrical Bose-Einstein Condensate 
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We present a calculation of a solitary wave propagating along a cylindrical Bose-Einstein trap, 
which is found to be a hybrid of a one-dimensional (ID) soliton and a three-dimensional (3D) vortex 
ring. The calculated energy-momentum dispersion exhibits characteristics similar to those of a mode 
proposed sometime ago by Lieb within a ID model, as well as some rotonlike features. 
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PACS numbers: 05.45.Yv, 03.75.Fi, 05.30.Jp 

Approximately forty years ago, Lieb and Liniger Q] 
employed the Bethe Ansatz to obtain an exact diagonal- 
ization of the model Hamiltonian that describes a one- 
dimensional (ID) Bose gas interacting via a repulsive 5- 
function potential. Based on this solution Lieb Q pro- 
posed an intriguing dual interpretation of the spectrum 
of elementary excitations, either in terms of the famil- 
iar Bogoliubov mode, or in terms of a certain kind of a 
particle- hole excitation which will hereafter be referred to 
as the Lieb mode. The energy-momentum dispersions of 
the Bogoliubov and the Lieb modes coincide at low mo- 
menta, where they both display the linear dependence 
characteristic of sound-wave propagation in an interact- 
ing Bose-Einstein Condensate (BEC), but the Lieb dis- 
persion is significantly different at finite momenta where 
it exhibits a surprising 2tt periodicity. Yet this new mode 
remained somewhat of a theoretical curiosity because of 
the absence of a physical realization of a strictly ID Bose 
gas. 

In a curious turn of events, the above subject re- 
emerged in connection with the solitary waves calculated 
analytically within a ID classical Gross-Pitaevskii (GP) 
model ||, f|], which were later shown to provide an ap- 
proximate but fairly accurate description of the quantum 
Lieb mode for practically all coupling strengths H |). 
More importantly, the same solitary waves motivated 
experimental observation of similar coherent structures 
in BEC of alkali-metal atoms through phase imprinting 
or phase engineering [0, Therefore, an opportunity 
presents itself for experimental realization of the Lieb 
mode. 

It is clear that a theoretical investigation based on ef- 
fective ID models often used to describe cylindrical traps 
will lead to the prediction of a Lieb mode, by arguments 
similar to those employed in the original ID classical 
model [|], ||. However, one should also question the sta- 
bility of the corresponding solitons within the proper 3D 
environment of a realistic trap ^, [l0|. For example, dark 
solitons created within a spherical trap were observed to 
decay into vortex rings plf . This observation was also 
supported by a numerical calculation in which an initial 
dark-soliton configuration is let to evolve in time accord- 
ing to the GP equation. 

Therefore, although the production of solitary waves 



through phase imprinting clearly suggests some distinct 
ID characteristics, it should be expected that a proper 
understanding of such waves must also account for 3D 
effects that are inevitably present in a realistic BEC. One 
could envisage a picture in which the actual solitary wave 
is a hybrid of a ID soliton and a 3D vortex ring. It is the 
aim of the present paper to make the above claim precise 
by calculating solitons propagating along a cylindrical 
trap without making a priori assumptions about their 
effective dimensionality. Our approach was motivated by 
the calculation of vortex rings in a homogeneous BEC 
due to Jones and Roberts ]12|] and a similar calculation 
of semitopological solitons in planar ferromagnets [ [l3| . 

Explicit results will be presented for a set of pa- 
rameters that roughly correspond to the experiment of 
Ref. m. Thus we consider a cigar-shaped trap filled 
with 8 'Tib atoms, with a transverse confinement fre- 
quency uo±_ = 2tt x 425 Hz and a corresponding oscilla- 
tor length a± = ^/h/muj± « 0.5/im. The coupling con- 
stant is written as U Q = A.iifi 2 a/m where a « 50 A is the 
scattering length. We make the approximation of an 
infinitely-long cylindrical trap with average linear density 
^ = 0.2 atoms/A and introduce the dimensionless combi- 
nations of parameters 7 = va = 10 and 7^ = va±_ = 10 3 . 
Finally, rationalized units are defined through the rescal- 
ings t— r^a^r, and ^ — > yjv ^ /a± . 

The energy functional extended to include a chemical 
potential is then given by 
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V** V* + p 2 + *#) 2 - 2/i dV, 

(1) 

where g = 4-7T7 and p 2 = x 2 + y 2 , and yields energy in 
units of jj_huj_\_ while the chemical potential \i is still 
measured in units of fiuo±_ . Similarly, the conserved linear 
momentum along the z axis given by the usual definition 
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is measured in units of his = j_\_(h/a±). In the second step 
of Eq. (||) we employed the usual hydrodynamic variables 



defined from \£ = 



An important first step in out calculation was to obtain 
accurate information about the ground state and the cor- 
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FIG. 1: Radial dependence of the solitary wave function at 
v = c/2 for various positive values of £ in steps of = 0.1. The 
result for negative £ is obtained through the symmetry rela- 
tions Re^(p,£) = Re^(p, -£) and Im^(p, £) = -Im^(p, -£). 
Distances are measured in units of aj_. 



responding linear (Bogoliubov) modes. The ground-state 
wave function ^ = ^o(p) was calculated as a stationary 
point of the energy functional (|l|) by a variant of a re- 
laxation algorithm jl4| and is normalized according to 
J °° 2i\pdp |^o | 2 = 1 to conform with our choice of ratio- 
nalized units. The chemical potential was found to be 
li = 6.4324 for 7 = 10. We have also recalculated the low- 
est branch of the Bogoliubov spectrum from which we 
extracted the speed of sound c = 1.77 in units of a_\_ujj_. 
This numerical estimate at 7 = 10 is consistent to within 
1% with the Thomas- Fermi approximation c = t^ 4 which 
was previously derived in a number of papers |y| [l6[ [l7| 
and was employed for the analysis of experimental data 

@. 

Thus we turn to the calculation of axially-symmetric 
solitary waves described by a wave function of the form 
\I/ = \I/(p, £) where (, = z — vt and v is the constant veloc- 
ity along the z axis. Such a wave function satisfies the 
stationary differential equation 
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±A* + ±p 2 * + g(***)*-iJL*, (3) 
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which is supplemented by the boundary conditions that 
\£ vanish for p — > oo, and — > |\£o(p)| f° r € ~> ±oo. 
The phase of the wave function is not fixed a priori at 
spatial infinity, except for a mild restriction implied by 
the Neumann boundary condition d^/dt; = imposed 
at £ — > ±oo mainly for numerical purposes. Finally, a 
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FIG. 2: Radial dependence of the local particle density 
n= |^| 2 for the four special cases discussed in the text. Con- 
ventions are the same as in Fig. [l]. 



solution of Eq. (||) must satisfy the virial relation 
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dV, 
(4) 

obtained by a standard scaling argument. 

Solutions of Eq. (Q) were obtained by a sophisticated 
Newton-Raphson iterative algorithm |l|, [l3| which will 
not be described here in any detail. For values of the 
velocity v near the speed of sound c, the calculated soli- 
tary wave is a rarefaction pulse that may be thought of 
as a mild soundlike disturbance of the ground state. The 
dominant features of the solitary wave are pronounced as 
the velocity is decreased to lower values and become rea- 
sonably apparent at v = c/2, as shown in Fig. [l] which pro- 
vides a complete illustration of the wave function through 
its real and imaginary parts. A partial but more trans- 
parent illustration is given in frame I of Fig. [| where we 
depict the radial dependence of the local particle density 
n= \^\ 2 for various values of £. It is clear that the density 
near the center of the soliton (£ = 0) vanishes on a ring 
with a relatively large radius R~2.8, thus a vortex ring 
is beginning to emerge. The features of the vortex ring 
become completely apparent, and its radius is tightened, 
as we proceed to yet smaller values of the velocity. A 
notable special case is the static (v = 0) vortex ring with 
radius R~ 1.8 illustrated in frame II of Fig. ||, which is 
far from being a completely dark (black) solitary wave. 

One would think that pushing the velocity v to nega- 
tive values would retrace the calculated sequence of soli- 
tary waves backwards. In fact, our algorithm continues 
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to converge to vortex rings of smaller radii until a criti- 
cal velocity v = — v — —0.475 c is encountered where the 
vortex ring achieves its minimum radius R m in — 0.8 and 
ceases to exist for smaller values of v. The terminal state 
at v = — vo is illustrated in frame III of Fig. ||. We have 
thus described a sequence of solitary waves with veloci- 
ties in the range — vo < v < c, which does not contain a 
black soliton. An equivalent sequence is obtained in the 
range — c < v < Vq simply by reversing the relative sign be- 
tween the real and imaginary part of the wave function, 
as is evident from Eq. (||). 

It is now important to calculate the energy-momentum 
dispersion. The excitation energy is defined as E = 
W — Wo where both W and Wq are calculated from 
Eq. ([!]) applied for the solitary wave \P(p, £) and the 
ground state ^o(p), respectively. The presence of the 
chemical potential in Eq. ([[]) provides the compensation 
that is necessary in order to compare energies of states 
with the same number of particles. Similarly, the rele- 
vant "physical" momentum is not the linear momentum 
V of Eq. (^|) but the "impulse" Q defined in a manner 
analogous to the case of a homogeneous gas || |l2| : 

Q = J \n-n )^dV = V-50 (5) 
2npdp no(p)[(l)(p,z = oo) - <j>(p,z = -oo)} , 



where no = |^o(p)| 2 is the ground-state particle density 
and 5(j) is the weighted average of the phase difference 
between the two ends of the trap. Here we simply postu- 
late the validity of definition (|[) and note that the corre- 
sponding group- velocity relation v = dE/dQ was satisfied 
to an excellent accuracy in our numerical calculation. On 
the other hand, the virial relation (||) was verified using 
the standard linear momentum V of Eq. (|^), as expected. 

The calculated dispersion E = E(Q) is depicted in 
Fig. H by a solid line along which we place the symbols 
I, II and III that correspond to the special cases of the 
solitary wave discussed above. At low Q the dispersion 
is linear, E = c |Q|, a feature that it shares with the Bo- 
goliubov dispersion. However, the energy now reaches a 
maximum at point II where v = and Q is slightly greater 
that 7r. Interestingly, the group velocity becomes nega- 
tive in the region (II, III) or, equivalently, the impulse is 
opposite to the group velocity. This rotonlike behavior 
could develop into a full-scale roton if the terminal point 
III actually turns out to be an inflection point beyond 
which the group velocity may begin to rise again. We 
have not been able to somehow continue our sequence 
of solitary waves beyond point III, nor have we ruled 
out such a possibility. In any case, the picture just de- 
scribed, together with the fact that the calculated radius 
of the vortex ring is monotonically decreasing along the 
sequence (I, II, III) comes close to the Onsager-Feynman 
view of a roton as the ghost of a vanished vortex ring 

The branch of the dispersion shown by a solid line in 
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FIG. 3: Energy E in units of j±(hco±) versus impulse Q in 
units of j±(h/a±). The solid line corresponds to the main se- 
quence of solitary waves discussed in the text, and the dashed 
line to the auxiliary sequence that contains a black soliton 
(point IV). 



Fig. H was calculated on the basis of the sequence of soli- 
tary waves with velocities in the range — vo <v<c. The 
branch for — c < v < vo may be thought to correspond 
to negative Q, as usual. However, there is an implicit 
2tt periodicity because of the appearance of the angle (j) 
in the definition of the impulse in Eq. (|5|). Therefore, 
a more natural representation of the spectrum seems to 
be obtained by restricting Q to the "first Brillouin zone" 
and thus completing Fig. || by its mirror image around 

Q = 7T. 

We have also carried out a detailed numerical cal- 
culation of the Bogoliubov dispersion for comparison. 
But the main point can be made here by assuming the 
model dispersion uj = \q\(c 2 + <7 2 /4) 1 / 2 where the fre- 
quency uo is measured in units of co± and the wave num- 
ber q in units of l/a±. Translated into the units em- 
ployed in Fig. pi, the model Bogoliubov dispersion reads 

E = \Q\(c 2 + 7lQ 2 / 4 ) 1/2 where a hu g e factor 7± ~ !0 6 
makes its appearance. Therefore, although the Bogoli- 
ubov and Lieb dispersions coincide at low Q, they diverge 
quickly at the scale of Fig. 0, a notable difference from the 
homogeneous ID model In other words, the two 

modes operate at rather different energy and momentum 
scales in a cylindrical trap. 

The last question we address is whether or not there 
exists an independent sequence of solitary waves that re- 
duces to a black soliton at v = 0. We used a model black- 
soliton configuration as input in our iterative algorithm 
which converged to the true static soliton illustrated in 
frame IV of Fig. ^| that is indeed a black soliton. We 
then incremented the velocity to either positive or neg- 
ative values with symmetrical results. For definiteness, 
we follow the solution to negative v and find that it ex- 
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FIG. 4: Contour levels of the local particle density n in a 
plane that contains the z axis and cuts across the cylindrical 
trap. The complete 3D picture may be envisaged by simple 
revolution around the z axis. Regions with high particle den- 
sity are bright while regions with zero density are black. The 
four special cases considered are the same as in Fig. 



ists only until one encounters the same critical velocity 
—vq discussed earlier in the text. The corresponding ter- 
minal state is precisely the same with the one reached 
through our original sequence of solitary waves and illus- 
trated in frame III of Fig. ^. For intermediate values of 
v the solution is again a vortex ring with constant radius 
at R = Rmin ~ 0.8. Accordingly, the calculated energy- 
momentum dispersion shown by a dashed line in Fig. || 
joins the original sequence through a cusp at the terminal 
point III, a situation that is reminiscent of the calculation 
of Jones and Roberts [Q. The same authors together 
with Putterman later argued j20| that the solitary waves 
that correspond to the upper branch are actually unsta- 
ble, a conclusion that might be valid in the present case 
as well. 

A summary view of our results is given in Fig. |. We 
have thus constructed a family of solitary waves which 
exhibit some quasi- ID features, such as the appearance 
of a nontrivial phase difference 5<fi that is important for 
phase engineering ||], but are otherwise bonafide 3D 
vortex rings. 
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